


eshock = -.01;

% Bivariate
irf2 = zeros(Tirf,3,varseps2);    %Each page corresponds to different model.
%
for i= 1:varseps2
    pvec = zeros(1,lags);         %treat as demeaned.
    xvec2nd = zeros(1,2*lags);
    if uselagflows ==1 
        %xvec = [4 lags of U/L and 4 lags of another outcome]
        xvec1st = zeros(1,2*lags);
    else
        %xvec = [4 lags of U/L]
        xvec1st = zeros(1,lags);
    end
    epsvec = zeros(1,lags);
    epsvec(1) = eshock;
    psum  = 0;
    %
    for t = 1:Tirf 
        pt= [  pvec, xvec1st, 0]*pcoef2(:,i) + (t==1)*eshock;
        ut= [epsvec, xvec2nd, 0]*ucoef2(:,i);
        yt= [epsvec, xvec2nd, 0]*ycoef2(:,i);
        xt= [ut;yt];
        %
        pvec_ = zeros(1,lags);
        pvec_(2:lags) = pvec(1:lags-1);
        pvec_(1) = pt;
        pvec = pvec_;
        %
        xmat = zeros(lags, 2);
        for j = 1:2
            xmat(2:lags,j) = xvec2nd(1+(j-1)*lags : j*lags-1);  % lags=4 and j=1 -> 1:3; j=2 -> 5:7 
            xmat(1,j) = xt(j);
        end
        xvec2nd = reshape(xmat,2*lags,1)';
        if uselagflows == 1
            xvec1st = xvec2nd;
        else
            xvec1st = xvec2nd(1:lags);
        end
        %
        epsvec_ = zeros(1,lags);
        epsvec_(2:lags) = epsvec(1:lags-1);
        epsvec_(1) = 0;
        epsvec = epsvec_;
        %
        if firstdiff ==1 
            % Cumulate responses of APL
            pt = pt + psum;
            psum = pt;
        end
        %
        irf2(t,1,i) = pt;
        irf2(t,2,i) = ut;
        irf2(t,3,i) = yt;
    end
end
%
%
% Trivariate
irf3 = zeros(Tirf,4,varseps3);      %Each page corresponds to different model.
%
for i= 1:varseps3
    pvec = zeros(1,lags);           %treat as demeaned.
    xvec2nd = zeros(1,3*lags);
    if uselagflows ==1 
        % xvec = [lags of U/L and f and one other outcome]
        xvec1st = zeros(1,3*lags); 
    else
        % xvec = [lags of U/L and f]
        xvec1st = zeros(1,2*lags);
    end
    epsvec = zeros(1,lags);
    epsvec(1) = eshock;
    psum  = 0;
    %
    for t = 1:Tirf 
        pt= [  pvec, xvec1st, 0]*pcoef3(:,i) + (t==1)*eshock;
        ut= [epsvec, xvec2nd, 0]*ucoef3(:,i);
        ft= [epsvec, xvec2nd, 0]*fcoef3(:,i);
        yt= [epsvec, xvec2nd, 0]*ycoef3(:,i);
        xt= [ut;ft;yt];
        %
        pvec_ = zeros(1,lags);
        pvec_(2:lags) = pvec(1:lags-1);
        pvec_(1) = pt;
        pvec = pvec_;
        %
        xmat = zeros(lags, 3);
        for j = 1:3
            xmat(2:lags,j) = xvec2nd(1+(j-1)*lags : j*lags-1);  % lags=4 and j=1 -> 1:3; j=2 -> 5:7 
            xmat(1,j) = xt(j);
        end
        xvec2nd = reshape(xmat,3*lags,1)';
        if uselagflows == 1
            xvec1st = xvec2nd;
        else
            xvec1st = xvec2nd(1:2*lags);
        end
        %
        epsvec_ = zeros(1,lags);
        epsvec_(2:lags) = epsvec(1:lags-1);
        epsvec_(1) = 0;
        epsvec = epsvec_;
        %
        if firstdiff ==1 
            % Cumulate responses of APL
            pt = pt + psum;
            psum = pt;
        end
        %
        irf3(t,1,i) = pt;
        irf3(t,2,i) = ut;
        irf3(t,3,i) = ft;
        irf3(t,4,i) = yt;
    end
end






